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ABSTRACT 


This  work  is  intended  to  provide  automated  methodology  for  processing  seismic  and  infrasound  data  from 
seismo-acoustic  arrays  and  apply  the  methodology  to  regional  networks  for  validation  with  ground  truth 
information.  In  the  initial  year  of  the  project  automated  techniques  for  detecting,  associating  and  locating 
infrasound  signals  were  developed.  Recently,  the  location  procedure  has  been  cast  into  a  Bayesian 
framework. 

Work  this  year  has  focused  on  assessment  of  sets  of  robust  regional  infrasound  signals  in  order  to  optimize 
the  location  algorithms.  Quantification  of  the  importance  of  model  and  empirical  travel  times,  backazimuth 
estimation  (including  errors)  and  signal  detectability  have  all  been  undertaken  and  are  reported.  Two 
specific  studies  are  reported. 

Infrasound  from  a  static  rocket  test  conducted  at  Promontory,  Utah  onl  1/01/2007,  was  observed  on 
USArray  seismographs  up  to  a  distance  of  600  km.  Signal  durations  were  broadly  consistent  with  the 
duration  of  the  test,  123  s,  and  signal  envelopes  correlated  with  variations  in  thrust  in  some  cases.  The 
effect  of  adding  Gaussian  noise  to  wind  profiles  is  investigated  for  the  purpose  of  phase  identification.  It  is 
shown  that  all  observations  can  be  predicted  through  the  combination  of  simulations  using  different  noise 
levels,  suggesting  that  a  probabilistic  approach  to  phase  identification  is  possible.  These  results  have 
important  implications  for  how  atmospheric  models  might  be  utilized  in  source  location  in  particular  how 
they  might  be  implemented  into  the  Bayesian  location  framework  of  Modrak  et  al.  (2010). 

The  deterministic  dependence  of  infrasound  propagation  on  the  season  and  path  environment  was 
quantified  by  the  analysis  of  more  than  1,000  repetitive  infrasonic  ground-truth  events  at  an  active,  open-pit 
mine  in  Korea  over  2  years.  Blast-associated  infrasonic  signals  were  analyzed  from  two  infrasound  arrays 
(CHNAR  and  ULDAR)  located  at  similar  distances  of  181  and  169  km,  respectively  from  the  source  but  in 
different  azimuthal  directions  and  with  different  path  environments.  The  CHNAR  array  is  located  to  the 
NW  of  the  source  area  with  primarily  a  continental  path,  whereas  the  ULDAR  is  located  E  of  the  source 
with  a  path  dominated  by  open  ocean.  As  a  result,  CHNAR  observations  were  dominated  by  stratospheric 
phases  with  characteristic  celerities  (group  velocities)  of  260  to  289  m/s  and  large  seasonal  variations  in  the 
travel  time,  whereas  data  from  ULDAR  consisted  primarily  of  tropospheric  phases  with  larger  celerities 
from  322  to  361  m/s  and  larger  daily  than  seasonal  variation  in  the  travel  time.  Ray  tracing  using 
atmospheric  models  incorporating  daily  weather  balloon  data  that  characterizes  the  shallow  atmosphere  for 
the  2  years  of  the  study  verifies  the  interpretation  of  these  observations.  Finally,  experimental  celerity 
models  that  included  seasonal  path  effects  were  constructed  from  the  long-term  dataset.  These  experimental 
celerity  models  were  used  to  constrain  travel  time  variations  in  infrasonic  location  algorithms  providing 
improved  location  estimates  as  illustrated  with  the  empirical  dataset. 
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OBJECTIVES 

There  are  two  primary  objectives  to  our  efforts  in  the  second  and  third  year  of  this  contract.  Each  objective 
builds  on  the  automatic  algorithms  for  detection,  association,  and  location  of  infrasound  events  developed 
in  the  first  year  of  this  project,  the  further  development  of  these  techniques  in  year  two  and  the  preliminary 
analysis  of  ground  truth  data  in  assessing  atmospheric  propagation  models. 

•  Validation  and  utilization  of  atmospheric  models.  Initial  efforts  in  this  contract  focused  exclusively  on 
infrasound  detection,  association,  and  location  without  requiring  the  use  of  atmospheric  models.  R&D 
this  year  has  focused  on  how  best  to  use  prediction  models  (comprising  the  combination  of  atmospheric 
and  propagation  models)  for  enhancing  source  localization.  In  this  paper  we  outline  two  different 
approaches:  a  stochastic  approach  and  a  deterministic  approach,  for  using  atmospheric  specifications  in 
event  location. 


•  Extended  analysis  of  infrasound  propagation  effects  using  ground-truth  data  from  Korea.  A 

primary  focus  in  the  last  year  has  been  in  the  careful  analysis  and  modeling  of  infrasound  observations 
associated  with  the  ground  truth  data  acquired  in  the  earlier  stages  of  the  study.  These  data  have  been 
compared  to  detailed  propagation  path  models  for  the  Korean  peninsula  and  provide  a  basis  for  assessing 
the  types  of  empirical  data  that  can  improve  location  estimates  based  exclusively  on  infrasonic  data  and 
seismic  and  infrasonic  data  in  combination. 


RESEARCH  ACCOMPLISHED 


A  Strategy  for  Enhancing  the  Bayesian  Infrasonic  Source  Location  (BISL)  Method 


The  BISL  framework,  outlined  in  Modrak  et  al.  (2010),  provides  a  robust  framework  for  infrasound 
location  without  requiring  a  prediction  model.  This  is  accomplished  in  part  by  defining  a  uniform  prior 
PDF  on  celerity  (group  velocity)  as  follows: 


1 16.67  0.22  km/s<  v  <  0.34  km/s 
[  0  otherwise 


This  Bayesian  prior  term  effectively  takes  a  worst-case  scenario:  we  have  no  expectation  of  the  celerity  for 
a  given  observed  arrival,  and  can  therefore  only  place  bounds  on  celerity  (noting  that  these  bounds  might 
be  more  restrictive  if,  for  example,  we  can  rule  out  a  thermospheric  return  on  some  basis).  In  the  absence  of 
a  reliable  prediction  model  that  is  capable  of  predicting  infrasound  phases,  this  appears  to  be  the  best  we 
can  do  for  now.  But  is  it  really? 

Prediction  Models 

Real  time  state-of-the-art  atmospheric  models,  such  as  the  Ground-to-Space  (G2S)  model  are  thought  to 
capture  the  day-to-day  and  hourly  variability  in  temperatures  and  winds  at  horizontal  resolutions  of 
approximately  1x1°.  However,  they  do  not  attempt  to  resolve  fine-scale  atmospheric  structures  such  as 
sub-grid  scale  gravity  waves,  turbulence,  and  wind  gusts.  It  has  been  demonstrated  that  the  combination  of 
such  atmospheric  models  with  existing  infrasound  propagation  models  do  not  always  reliably  predict 
observations  (e.g.,  Arrowsmith  et  al.;  2007;  Evers,  2009).  This  is  a  particular  problem  in  the  zone-of- 
silence,  a  region  that  extends  from  approximately  50  to  200  km.  Two  factors  influence  this  limitation: 

(1)  the  fidelity  of  atmospheric  models,  and  (2)  the  fidelity  of  propagation  models. 

For  routine  infrasound  location,  we  require  a  good  blend  of  simplicity  and  physical  thoroughness,  to 
balance  the  need  to  be  useful  in  operational  processing  and  capture  sufficient  propagation  physics.  As 
outlined  by  Drob  et  al.  (2010),  the  TauP  method  (e.g.,  Garces  et  al.,  1998;  Drob  et  al.,  2010)  provides  a 
good  balance  of  simplicity  and  physical  thoroughness  to  make  it  the  propagation  model  of  choice  for 
automatic  location  algorithms.  Other  methods  -  such  those  based  on  Parabolic  Equations  or  Normal  Modes 
(see  Whitaker  and  Norris,  2008,  for  a  summary)  -  account  for  additional  wave  physics,  such  as  diffraction 
effects,  but  are  impractical  for  current  automated  techniques  due  to  the  high  computational  costs. 

Using  Prediction  Models  in  BISL 

The  simple  question  is  this:  Given  some  smoothed  realization  of  the  atmosphere  at  a  given  place  and  time, 
how  can  we  use  that  information  to  better  constrain  the  group  velocity?  With  the  knowledge  that  we  cannot 
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reliably  predict  arrivals  deterministically,  we  need  to  somehow  account  for  the  unknown  effects  of  gravity 
waves,  turbulence,  and  wind  gusts  -  as  well  as  limitations  in  propagation  physics  -  in  a  stochastic  way.  A 
simple  way  to  do  this  is  to  add  noise  to  G2S  wind  profiles,  constrained  by  empirical  data.  We  note  here  that 
Gibson  et  al.  (2009)  have  used  spectral  models  for  gravity  waves  to  perturb  atmospheric  models  with  some 
success,  using  a  Time  Domain  Parabolic  Equation  method.  Since  our  strategy  is  to  account  for  multiple 
sources  of  error  for  event  location,  a  simple  Gaussian  noise  model  is  used  here.  Incorporating  additional 
constraints  such  as  the  gravity  wave  spectral  model  of  Gardner  (1995)  could  refine  this  model. 

Based  on  the  fact  that  the  dense  network  of  seismic  stations  in  the  USArray  can  be  used  for  detecting 
infrasound  (e.g.,  de  Groot-Hedlin  et  al.,  2008),  we  can  leverage  this  unprecedented  sampling  of  the 
infrasound  wavefield  to  test  the  effect  of  adding  different  levels  of  noise  to  the  zonal  and  meridional  wind 
components,  and  determine  what  standard  deviation  of  noise  is  needed  to  predict  the  observations.  A  static 
rocket  test,  conducted  at  Promontory,  Utah,  at  19:00  UT  on  1 1/01/2007  was  well  recorded  by  USArray  and 
easy  to  identify  on  the  single-channel  seismic  data  owing  to  its  unique  source-time  function  of  duration 
equal  to  123  s  (Figure  1). 


Azimuth=135 


Figure  1.  Record  section  taken  along  a  great-circle  azimuth  of  135°.  All  stations  within  20°  of  the 

great-circle  azimuth  are  projected  onto  the  great-circle  path.  Vertical  component  waveforms 
are  plotted,  bandpass  filtered  from  3-10  Hz.  Red  lines  correspond  to  group  velocities  of  0.34 
km/s,  0.29  km/s,  and  0.22  km/s. 

We  add  noise  to  the  wind  components  only,  since  the  uncertainties  in  temperature  in  the  G2S  model  are 
lower  (Drob  et  al.,  2010).  The  implementation  is  simple:  Mathematically,  u(z)  =  w0(z )  +  7J  and 
v(z)  =  V0(z)  +  //,  where  U0  and  V0  are  the  unperturbed  G2S  wind  profiles  and  //  is  the  noise  vector. 
Figure  2  illustrates  the  addition  of  noise  by  comparing  zonal  and  meridional  winds  with  and  without 
Gaussian  noise  with  a  standard  deviation  of  5  m/s. 

A  comparison  between  the  observations  and  predictions,  for  different  levels  of  added  noise,  is  shown  in 
Figure  3.  Gaussian  noise  with  a  equal  to  5  m/s,  10  m/s,  and  20  m/s  have  been  added  to  the  unperturbed 
wind  profiles  in  Figure  2.  The  uncertainties  used  are  within  the  approximate  range  of  estimated 
uncertainties  in  G2S  models  (Drob  et  al.,  2010),  although  the  variation  in  uncertainty  with  height  is  not 
accounted  for  here.  The  predictions  estimated  here  do  not  include  direct  arrivals  --  only  ducted  arrivals  are 
considered.  Furthermore,  multiple  bounces  are  not  considered  in  this  study,  eliminating  the  added 
complexity  of  accounting  for  topography  (Arrowsmith  et  al.,  2007).  Multiple  atmospheric  realizations  have 
been  computed  at  each  noise  level,  with  similar  ray-tracing  results  obtained  each  time.  The  results 
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presented  below  correspond  to  single  realizations,  but  these  are  broadly  representative  of  multiple  runs  in 
terms  of  the  highlighted  regions  where  predictions  match  (see  Figure  3). 
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Figure  2.  Profiles  of  the  ID  temperature  (left  panel)  and  winds  (right  panel)  extracted  from  the  G2S 
model  for  19:00  on  11/01/2007  at  the  test  location.  Unperturbed  zonal  and  meridional  winds 
are  represented  by  blue  lines,  while  perturbed  wind  profiles  (using  £7=5  m/s)  are  shown  in 
red. 

The  predicted  bounce-points  (locations  where  rays  are  predicted  to  hit  the  ground),  which  are  shown  in 
Figure  3,  assume  an  omni-directional  source  and  do  not  account  for  the  known  source  directivity.  The 
bounce-points  were  predicted  by  shooting  a  suite  of  rays  at  azimuths  of  0-355°,  and  launch  angles  of 
0-90°,  with  5°  spacing  in  both  azimuth  and  launch  angle.  With  no  added  noise,  the  geographic  distribution 
of  bounce-points  reflects  the  strong  positive  zonal  wind  in  the  stratosphere,  common  at  this  latitude  during 
the  winter  months,  which  enhances  the  stratospheric  duct  to  the  east  and  destroys  the  stratospheric  duct  to 
the  west  (noting  that  stratospheric  group  velocities  are  typically  0.28  km/s-0.31  km/s). 

The  results  shown  in  Figure  3  highlight  that,  although  some  arrivals  are  not  predicted  using  the  unperturbed 
G2S  profdes,  they  can  be  predicted  when  noise  is  added  to  the  profdes.  For  example,  observations  with  Vg 

=  0.28  km/s  to  the  southwest  of  the  test  are  only  predicted  with  o  =  5-10  m/s;  observations  with  Vg  <  0.27 

to  the  northeast  of  the  test  are  only  predicted  with  a  =  20  m/s.  As  would  be  expected,  with  increasing 
noise,  returns  are  predicted  at  ever-closer  distances.  Flowever,  due  to  the  strong  zonal  winds  that  favor 
propagation  to  the  east,  few  returns  are  predicted  towards  the  west  at  these  noise  levels,  with  the  exception 
of  thermospheric  returns. 

The  full  suite  of  observations  can  only  be  explained  through  the  combination  of  simulations  with  varying 
noise  levels.  This  detail  is  important  because  it  implies  that,  for  phase  identification  (and  thus  also  for 
location  purposes),  a  probabilistic  approach  is  necessary.  For  example,  an  observed  arrival  may  be 
explained  as  either  a  stratospheric,  or  tropospheric,  arrival  depending  on  the  standard  deviation  of  the  wind 
perturbation  addition. 


670 


2010  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


Our  contention  is  that  if  we  can  reliably  predict  arrivals  through  the  addition  of  noise  to  atmospheric  wind 
models,  the  simple  uniform  prior  of  Modrak  et  al.  (2010)  may  be  refined  using  a  stochastic  prediction 
model.  By  the  principle  of  reciprocity,  the  Tau-P  equations  can  be  calculated  from  the  receiver  perspective 
by  integrating  the  ray  equations  with  a  negative  time  step.  Thus,  by  tying  the  standard  deviation  of  wind 
perturbations  to  probability  (where  more  probable  predictions  are  based  on  smaller  standard  deviations  of 
additive  noise),  we  can  develop  a  non-uniform  prior  PDF,  for  a  given  signal  backazimuth  and  arrival  time, 
based  on  Tau-P  predictions.  In  effect,  we  can  directly  relate  celerity  to  the  standard  deviation  of  wind 
perturbations,  and  subsequently  relate  wind  noise  probabilities  to  group  velocity  probabilities.  Further 
research  is  needed  to  develop  this  R&D  concept,  implementing  prediction  model  simulations  into  BISL2.0. 
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Figure  3.  Plots  comparing  predictions  (open  circles)  and  observations  (filled  triangles)  for  different 
levels  of  noise  added  to  the  atmospheric  data.  Black  dots  represent  stations  that  did  not 
record  the  test.  The  predictions  and  observations  are  plotted  on  the  same  color  scale,  which 
represents  group  velocity.  Thermospheric  predictions  have  been  removed  for  clarity. 
Observations  enclosed  by  black  lines  are  uniquely  predicted  for  that  given  noise  level. 
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Experimental  Characterization  of  Seasonal  Variations  in  Infrasonic  Travel  Times  in  the  Korean 
Peninsula  with  Implications  for  Infrasound  Location 

The  purpose  of  this  component  of  our  work  is  to  experimentally  characterize  infrasonic  travel  time  and 
construct  quantitative  celerity  models  that  depend  on  the  season  and  propagation  path.  In  support  of  this 
deterministic  approach  to  travel-time  assessment  we  collected  nearly  daily  repetitive  infrasonic 
observations  over  a  two-year  time  period  from  ground-truth  events  at  an  active  mine.  These  events  included 
more  than  1,000  surface  blasts  that  occurred  daily,  with  multiple  events  on  active  days,  over  the  2  years. 
Variations  in  infrasound  propagation  based  on  the  data  are  quantified  at  the  two  infrasound  arrays  CHNAR 
and  ULDAR  in  Korea  (Figure  4),  both  located  at  similar  distances  from  the  mine  but  with  different 
propagation  environments  -  one  path  dominated  by  the  continent  and  the  other  by  the  ocean. 


Figure  4.  Locations  of  the  open-pit  mine  (circle),  two  infrasound  arrays  (stars)  and  weather  balloon 
station  (square)  used  in  this  study.  A  photograph  of  one  in-mine  seismo-acoustic  station  is 
shown  in  the  lower  right. 

Array  processing  based  on  the  progressive  multi-channel  correlation  (PMCC)  (Cansi,  1995),  is  used  to 
detect  and  analyze  the  blast-associated  infrasonic  signals  and  to  construct  a  celerity  model  for  each  path. 
Characteristics  of  the  two  sets  of  experimental  celerity  models  include  seasonal  dependent  path  effects  to 
each  array.  Ray  tracing  using  long-term  atmospheric  models  supplemented  by  shallow  atmosphere 
constraints  provided  by  weather  balloon  data  was  performed  for  comparison  with  the  observations. 

Variations  of  infrasound  propagation.  During  the  two  years  of  in-mine  sensor  operation,  1,066  blasts 
were  detonated,  more  than  one  blast  per  day  on  average.  97%  of  the  blasts  occurred  at  or  around  03  UTC 
(12  pm  local  time).  The  origin  time  and  approximate  location  of  each  event  was  estimated  using  the  seismic 
and  infrasound  observations  from  the  in-mine  stations.  Infrasonic  array  data  were  analyzed  using  the 
PMCC  method,  in  order  to  detect  blast-associated  infrasonic  signals.  The  wave  parameters  of  each  detected 
signal,  including  the  arrival  time  and  azimuth,  were  then  estimated  at  each  array.  Figure  5  displays  the 
waveforms  and  analysis  results  for  one  set  of  the  ground-truth  events.  The  left  figure  shows  near-field 
seismic  and  acoustic  signals  recorded  at  the  two  in-mine  stations  (ACE3  and  ACE7).  The  acoustic  signals 
in  this  case  were  clipped  due  to  instrument  dynamic  range  and  the  close  distance  to  the  10-ton  blast.  We 
simply  used  the  center  of  the  main  blasting  zone  as  a  source  location  for  all  blasts.  PMCC  results  are  also 
shown  in  the  figure.  WinPMCC  v2.1  was  used  to  make  the  regional  wavefield  estimates  focusing  on  a 
phase  velocity  range  of  0.25-0.40  km/s  using  a  window  length  20  s,  and  a  frequency  band  0.5-5  Hz.  After 
signal  association,  we  compiled  travel  times  and  wave  parameters  determined  from  PMCC. 
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Figure  5.  Example  data  for  a  blast  charged  with  10  tons  of  explosives,  (left)  The  four  waveforms  from 
the  blast  are  the  near-field  seismic  and  acoustic  signals  recorded  in  the  mine,  (right)  Boxes 
mark  the  blast-associated  infrasound  signals  recorded  at  the  two  regional  arrays  with  the 
upper  ten  waveforms  from  CHNAR,  excluding  one  bad  channel  (CHN32),  and  next  seven 
waveforms  from  ULDAR.  The  corresponding  PMCC  results  display  the  estimates  of 
azimuth  (upper)  and  apparent  velocity  (lower)  for  the  detected  signals.  Measured  travel 
times  for  both  infrasonic  signals  are  listed  in  the  lower  right. 


Figure  6  documents  the  variation  of  infrasound  travel  times  and  azimuths  to  the  two  arrays  for  all  the 
ground-truth  events  that  were  observed  over  the  2  years.  Despite  similar  propagation  distances,  the  average 
travel  times  were  significantly  different  to  the  two  arrays,  with  arrivals  at  ULDAR  on  average  150  seconds 
earlier  than  those  at  CHNAR.  For  CHNAR,  the  celerities  calculated  from  the  distances  and  travel  times 
ranged  from  260  to  289  m/s,  which  corresponds  to  the  characteristic  value  of  stratospheric  returns 
(Evers  and  Haak,  2007). 
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Month  of  year,  2007-2009 

Figure  6.  Variations  of  travel  times  (circles)  and  azimuths  (triangles)  estimated  at  both  arrays  from 
the  ground-truth  events  recorded  over  the  2  years  of  the  study.  The  dotted  lines 
superimposed  on  the  observed  travel  times  are  sinusoidal  celerity  models  that  were  lit  to 
the  observations.  Note  that  ULDAR  was  installed  in  December  2007  and  so  there  is  a 
shorter  time  period  of  observation. 


In  addition,  the  travel  times  showed  systematic  seasonal  variations  with  signals  arriving  up  to  68  seconds 
faster  in  the  summer.  Thus,  the  celerity  showed  a  smoothly  varying  cyclical  change  with  time,  with 
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velocities  of -0.26,  -0.29,  and  -0.27  km/s  in  the  winter-spring,  summer,  and  fall  seasons,  respectively. 

This  systematic  variation  in  observed  celerity  is  consistent  with  the  contention  that  the  travel  times  of 
infrasound  signals  are  strongly  dependent  on  the  seasonal  ambient  sound  speed  along  the  path.  In  contrast, 
the  infrasound  signals  arriving  at  ULDAR  showed  relatively  large  celerities  (322  to  361  m/s),  characteristic 
of  tropospheric  returns.  The  travel  times  to  ULDAR  show  little  seasonal  variation  but  document  larger 
daily  variations  of  up  to  56  second  in  the  travel  times  between  days.  These  observations  can  be  partially 
attributed  to  the  dominant  effect  of  the  ocean  on  temperatures  at  the  ocean-atmosphere  boundary. 

Path  variations  also  impact  the  ability  to  detect  the  arrivals.  Signals  at  CHNAR  were  most  likely  to  be 
detected  in  the  summer.  Many  signals  were  not  detected  during  the  fall-winter-spring  seasons,  even  though 
strong  source  signals  were  identified  at  the  in-mine  stations  for  these  times.  Detection  at  ULDAR  was  low 
in  summer  (August)  but  stable  during  the  other  seasons.  Thus,  the  seasonal  detection  at  the  two  arrays  was 
complementary,  as  shown  in  Figure  6  as  result  of  the  azimuthal  differences  of  the  two  arrays  relative  to  the 
mine.  This  variation  in  detection  can  be  partially  explained  by  seasonal  changes  in  the  atmospheric  wave 
duct  attributable  to  known  seasonal  wind  in  and  around  the  Korean  Peninsula  with  winds  northwesterly  in 
the  winter  and  southeasterly  in  the  summer.  In  contrast  to  the  variation  of  travel  times,  the  azimuth 
estimates  showed  little  variation  with  season  at  either  array  (Figure  7).  The  mean  and  variance  of  the 
azimuth  estimates  were  1 15.6±1.6°  at  CHNAR  (1 15.8°  of  true  azimuth)  and  271.5±3.1°  at  ULDAR 
(271.7°),  respectively.  The  insensitivity  of  backazimuth  to  seasonal  variations  is  consistent  with  the 
relatively  weak  meridional  winds  (5.9%  and  27.2%  of  annual  average  zonal  wind  speed  at  10  and  30  km 
above  the  mine,  respectively,  estimated  from  Horizontal  Wind  Models  (Drob  et  al.,  2008)  in  the  region, 
which  are  normal  to  the  dominantly  EW  propagation  paths  in  this  study. 


Figure  7.  Distribution  of  backazimuth  estimates  at  CHNAR  (left)  and  ULDAR  (right)  for  all  the 
ground  truth  events  observed  at  the  two  arrays.  The  blue  line  represents  the  true 
backazimuth  to  the  mine  where  the  events  occurred. 

Modeling  of  Infrasonic  Propagation.  Atmospheric  ray  tracing  was  performed  for  all  the  observations 
over  the  2  years  of  the  study  using  a  combination  of  empirical  atmospheric  models  (Horizontal  Wind 
Model/Mass  Spectrometer,  Incoherent  Scatter:  HWM07  (Drob  et  al.,  2008)/MSISE-90  (Hedin  et  al.,  1991) 
for  the  upper  atmosphere  and  weather  balloon  data  for  the  lower  atmosphere.  The  lower  part  of  the 
effective  sound  velocity  profiles  (ESVP)  used  for  ray  tracing  (Hagerty  et  al.,  2002)  was  based  on  the 
weather  balloon  data  from  the  Sokcho  station,  which  is  closer  to  both  the  mine  and  arrays  than  other 
stations  operated  by  the  KMA.  The  measurement  height  of  the  balloon  data  depended  on  the  day,  but  data 
were  usually  available  to  an  approximate  altitude  of  30  km  and  included  temperature,  wind  direction,  and 
wind  speed.  Data  acquired  at  00  UTC  was  used  for  the  modeling  because  it  is  within  3  hours  of  the  average 
time  of  the  blasts.  Ray  tracings  at  a  0. 1°  intervals  were  computed  for  each  day  of  the  experimental  period 
and  included  models  for  726  days  at  CHNAR  and  484  days  at  ULDAR  (smaller  number  of  days  for 
ULDAR  reflects  its  more  recent  installation,  December  2007).  When  trajectories  of  a  ray  passed  within  1 
km  (capture  radius)  of  the  array,  the  ray  was  identified  as  an  eigenray  for  the  infrasonic  phase.  Of  the 
eigenrays  computed  on  a  given  day,  the  ray  with  highest  turning  height  was  selected  as  a  representative 
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eigenray  for  the  day.  In  case  of  CHNAR.  the  total  rays  from  the  simulations  were  grouped  into  four 
dominant  classes:  (1)  in  43.8%  of  the  cases  there  exists  a  thermospheric  phase  but  CHNAR  is  located  in  a 
acoustic  shadow  zone  for  the  rays,  It_shadow\  (2)  in  36.9%  of  the  cases  there  exists  a  stratospheric  phase 
but  CHNAR  is  again  located  in  a  shadow  zone,  Is_shadow\  (3)  in  16.7%  of  the  cases  there  was  a 
tropospheric  phase  that  was  predicted  to  return  to  CHNAR,  Iw_eigeir,  and  (4)  in  2.6%  of  the  cases  there 
was  predicted  a  stratospheric  phase  that  is  predictable  to  return  to  CHNAR,  Is_eigen  .  In  case  of  ULDAR, 
four  dominant  cases  were  identified  and  included  It_shadow  (9.9%),  Is_shadow  (5.4%),  Iw_shadow 
(3.1%),  and  Iw_eigen  (81.6%).  Unlike  the  other  shadowing  phase,  the  Iw_shadow  phases  jumped  ULDAR 
with  bounces  at  the  surface. 

Figure  8  displays  the  variation  of  ESVP  for  the  2  years  of  the  study  along  with  the  ray-tracing  results.  The 
refraction  altitudes  of  each  daily  eigenray  from  the  mine  to  the  individual  arrays  are  superimposed  on  the 
ESVP. 


Month  of  year.  2007-2009  Month  of  year.  2007-2009 

Figure  8.  Ray  tracing  results  using  effective  sound  velocity  profiles  (ESVP)  for  the  path  from  the 
mine  to  CHNAR  (left)  and  to  ULDAR  (right)  for  the  2  years  of  the  study.  White  dots 
indicate  boundary  heights  between  balloon  data  and  the  empirical  data  used  for  the  ray 
tracing.  Black  dots  indicate  the  lowest  possible  refraction  heights  where  the  sound  speed  is 
larger  than  that  at  source.  White  circles  indicate  turning  points  in  the  atmosphere  from  the 
predicted  eigenrays  for  each  day  at  00  UTC. 

For  CHNAR,  ray  tracing  predicted  a  smaller  number  of  eigenrays  representative  of  stratospheric  returns 
than  that  documented  in  Figure  6.  Although  many  stratospheric  returns  were  predicted  from  ray  tracing 
(the  black  dots  around  40  km  altitude),  CHNAR  was  in  the  so-called  zone  of  silence  in  most  of  these 
prediction.  This  comparison  implies  that  the  ESVP,  particularly  the  empirical  models  with  stratospheric 
returns,  were  not  sufficient  to  predict  the  real  observations  or  that  the  capture  radius  was  too  small,  even 
though  the  predictions  for  the  Is_eigen  are  increased  by  widening  the  capture  radius,  for  example  4.7%  for 
5  km  and  14.5%  for  10  km  of  capture  radius.  Most  Is_eigen  phases  at  CHNAR  were  predicted  during  the 
winter  (see  January  of  2008  and  2009)  with  refraction  altitudes  corresponding  to  heights  sampled  by  the 
balloon  data.  These  predictions  can  be  explained  by  not  only  the  slow  sound  speed  at  the  ground  but  also 
the  slightly  higher  velocity  channel  at  stratospheric  heights  during  this  season.  The  higher  level  of  signal 
detection  predicted  for  the  mid- winter  (January)  in  comparison  to  the  spring  and  fall  seasons  was  consistent 
with  the  observations  detailed  in  Figure  6.  The  calculations  also  predicted  many  Iw_eigen  phases  that 
refracted  at  very  low  altitudes  and  experienced  bounces  at  the  ground  surface  but  the  corresponding  phases 
were  not  clearly  observed  in  the  ground  truth  data.  This  second  result  suggests  that  the  Iw  phase  may  have 
limited  horizontal  propagation  along  the  continental  path  which  includes  significant  topographic  variation 
between  the  mine  and  CHNAR.  In  contrast,  many  eigenrays  returned  to  ULDAR  by  refraction  at  altitudes 
lower  than  20  km.  These  altitudes  were  in  the  range  where  the  balloon  data  were  collected. 

Infrasound  Location  of  Ground  Truth  Events.  Infrasonic  signals  from  130  blasts  were  detected  by  both 
arrays.  In  this  study,  we  calculated  two  different  location  sets  for  these  selected  blasts.  The  first  set  of 
location  estimates  were  determined  by  assuming  common  stratospheric  phase  for  both  arrays  with  constant 
celerity  of  0.3  km/s  for  all  seasons  in  the  least-square  method.  The  second  set  of  location  estimates  were 
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calculated  using  simple  sinusoidal  celerity  model  which  were  fit  to  the  observation  (dotted  lines  in 
Figure  6),  and  includes  seasonal  dependences  as  described.  The  location  results  are  summarized  in  Figure  9 
with  the  blue  squares  representing  the  first  set  of  locations  and  the  red  ones  the  second  location  set  using 
the  seasonal  dependent  celerity  models.  95%  confidence  ellipses  for  each  location  set  are  also  included  in 
the  figure.  The  centroid  for  the  location  of  the  first  set  is  biased  about  17.8  km  to  the  east  of  the  mine  due  to 
the  discrepancies  between  the  assumed  constant  celerity  and  observed  celerities  shown  in  Figure  6.  The 
average  location  estimates  using  the  experimentally  determined  celerity  models  that  included  seasonal  path 
effects  showed  only  a  1.7  km  offset  from  the  mine  location.  This  result  suggests  that,  for  infrasonic  source 
location  at  regional  distances,  the  dependence  of  infrasound  propagations  on  seasons  and  path 
environments  should  be  quantified  and  its  effects  incorporated  into  the  infrasound  location  method  in  order 
to  produce  improved  location  estimates.  In  the  figure,  black  circles  indicate  seismic  locations  for  the  blasts 
extracted  from  the  seismological  catalog  of  the  KEMS.  Of  the  total  ground  truth  dataset  used  in  this  study, 
the  KEMS  detected  589  events  (open  circles)  with  local  magnitude  ranging  from  0.5  to  1.8  ML.  Yellow 
circles  in  the  figure  indicate  the  100  seismic  locations  which  correspond  to  the  130  blasts  located  using 
only  the  infrasound  signals  as  described  above.  It  is  interesting  to  note  that  there  were  30  events  in  the 
infrasound  set  that  did  not  produce  seismic  events  in  the  database.  Comparing  the  infrasonic  location 
with  the  seismic  location,  suggests  that  the  infrasonic  location  including  seasonal  path  effects  are  at 
least  comparable  to  seismic  location  in  case  of  small-magnitude  surface  explosions  at  regional 
distances. 


Figure  9.  Two  infrasonic  locations  results  for  selected  ground  truth  events  using  constant  celerity 

model  (blue  squares)  and  sinusoidal  model  (red  squares).  Circles  indicate  seismic  locations 
of  the  ground  truth  events  determined  by  the  KEMS. 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  stochastic  study  suggests  that  one  can  reliably  predict  arrivals  through  the  addition  of  noise  to 
atmospheric  wind  models.  The  implication  is  that  the  simple  uniform  prior  of  Modrak  et  al.  (2010)  may  be 
refined  using  a  stochastic  prediction  model.  By  the  principle  of  reciprocity,  the  Tau-P  equations  can  be 
calculated  from  the  receiver  perspective  by  integrating  the  ray  equations  with  a  negative  time  step.  Thus,  by 
tying  the  standard  deviation  of  wind  perturbations  to  probability  (where  more  probable  predictions  are 
based  on  smaller  standard  deviations  of  additive  noise),  we  can  develop  a  non-uniform  prior  PDF,  for  a 
given  signal  backazimuth  and  arrival  time,  based  on  Tau-P  predictions.  In  effect,  we  can  directly  relate 
celerity  to  the  standard  deviation  of  wind  perturbations,  and  subsequently  relate  wind  noise  probabilities  to 
group  velocity  probabilities. 

The  complementary  two-year  deterministic  study  of  infrasound  wave  propagation  has  utilized  ground  truth 
events  collected  at  an  active  mine  in  Korea.  The  dataset  included  1,066  blasts  and  highlights  the  variation 
of  infrasonic  propagation  due  to  the  seasonal  variation  of  atmospheric  structures  and  propagation  path 
environments.  The  blast-associated  infrasonic  signals  were  analyzed  at  two  distant  infrasound  arrays,  which 
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were  located  at  similar  distances  (zone  of  silence)  from  the  mine  but  pointed  in  different  azimuth  directions 
over  different  path  environments  (continental  and  open  ocean).  Infrasonic  waves  traveling  NW  toward 
CHNAR  over  the  continent  propagated  as  guided  waves  between  the  ground  and  stratosphere,  with  a 
characteristic  celerity  range  of  260-289  m/s.  These  observations  documented  seasonal  cyclical  variations  in 
travel  time.  The  infrasonic  waves  traveling  E  to  ULDAR  also  propagated  as  guided  waves,  but  observed 
celerities  indicated  that  wave  ducts  were  formed  at  relatively  low  heights  in  the  troposphere.  The  celerity 
range  was  correspondingly  large  (322-361  m/s),  even  though  the  propagation  distances  to  both  arrays  were 
similar.  ULDAR  showed  more  daily  variation  in  travel  time  than  seasonal  variation.  The  lower  refractions 
were  compared  and  verified  by  ray  tracings  with  balloon  data  constraining  the  lower  atmosphere.  In 
particular,  the  observations  at  ULDAR  support  the  existence  of  an  ephemeral  ‘SOFAR’  layer  in  the 
atmosphere  (Herrin  et  al.,  2006)  throughout  the  year  over  the  open  ocean.  Finally,  the  dataset  illustrates  that 
a  celerity  model  that  captures  the  seasonal  path  effects  can  substantially  improve  infrasound  location 
procedures.  The  infrasound  only  location  resulted  in  improved  locations  when  compared  routine  seismic 
locations  for  these  small-magnitude  surface  explosions  observed  at  regional  distances. 
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